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ABSTRACT 


In this study, the time temperature profile of a missile exposed to fire in a 
compartment adjacent to the missile magazine is examined. The study required the 
development of a heat transfer model based on the geometry and thermophysical 
properties of a new concept for a vertical launching system, the Concentric Canister 
Launcher (CCL). Different fire scenarios are analyzed by the model to predict the time it 
takes to reach a critical value or “cook-ofF’ temperature of the missile’s propellant and 
explosives. 
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I. INTRODUCTION 


A. BACKGROUND 

The history of naval combat has shown that there is a demonstrated need to 
understand how shipboard fires, and their containment, affect ship survivability. Credible 
fire threat models are particularly helpful in assessing the adequacy of fire protection 
systems, and in providing baseline data to make sound damage control decisions to 
combat conflagrations which may jeopardize munition magazines. Fires initiated by the 
penetration and detonation of Excocet missiles, and sustained with their unspent missile 
propellent, were experienced by the USS STARK (FFG 31) in the Arabian Gulf and 
British ships during the Falkland conflict. In all instances, ship’s survivability was 
seriously jeopardized. In the Falklands, a British ship was abandoned and later lost. In 
each instance, munitions magazines were threatened by fire that spread into adjacent 
compartments. Drastic measures were required to cool the munitions, and prevent the 
detonation of high explosive material and the ignition of solid propellent. [Ref 1] 

B. CCL CONCEPT 

Naval Surface Warfare Center Dahlgren Division (NSWC) [Ref 2] is currently 
examining a new concept for a future vertical launching system, the Concentric Canister 
Launcher (CCL). This new launcher will provide greater firepower, flexibility, and lower 
cost based on its modular design, than the current MK 41 VLS. The CCL was originally 
proposed to replace the MK 26 Guided Missile Launching System (GMLS) aboard CGN 
36, DDG 993, CG 47-51, and non VLS, DD 963 classes to enhance capability with 
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increased magazine capacity. Additionally, the CCL would be adaptable to forward fit 
new ships such as future DDG-Sl's, SC-21, Arsenal and Swath hull designs. 

A schematic of the CCL is given in Figure 1. It consists of two concentric 
cylinders joined by sets of dual longerons. One end is open and the other is sealed with a 
hemispherical cap. The annulus between the inner and outer cylinder provides gas 
management by ducting discharge exhaust combustion gases at launch. Each canister 
carries one missile and fits into a ship’s weapon module. The weapons module consists of 
a test tube type rack in which each missile canister is held by a shock collar at the main 
deck. Schematics of the ship’s weapon module and shock collar are given in Figures 2 
and 3, respectively. The potential advantages of this concept over the current MK 41 are: 

• Round pressure vessels and one shot design allow for the lightest possible 
configuration. This frees up ship’s weight restrictions and allows for more 
missiles to be carried. 

• Missiles fire vertically, horizontally, or at any angle in between. Thus, 
allowing CCL’s to be outfitted in geometrically constraining spaces such as 
the narrow hull designs of the swath and surface-effects ships. 

• A shock collar mitigates the underwater shock at the main deck rather than 
the inner bottom or platform level. 

• A wide assortment of different weapons loadouts, such as ships self 
defense weapons (Seasparrow AAW missiles), electronic warfare devices 
(chaff and flares), torpedoes (MK 46 and 50), remote sensing devices 
(sonobuoys) in addition to the current capabilities of the MK 41 VLS 
(TEAM, SLAM, NATACMS, SM-2, etc ). 
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Figures. Shock Collar. 
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The CCL, although different from the MK 41 VLS, is not intended to replace it, but 
merely provides ship builders with a new set of options for launching a wide variety of 
weapons. 

A CCL module sized to fit in the after VLS space of a DDG-51 ship contains up to 
80 Tomahawk missiles. This weapons loadout contains approximately 24,000 lb. of solid 
propellant, 30,000 lb. of explosives, and 35,300 lb. of liquid fuel [Ref 1]. Thus the fires 
experienced by British ships in the Falklands conflict and the USS STARK (FFG 31) 
involving this much high energy material has provided a vital insight into the dangers and 
difficulties of fire engulfing the spaces surrounding a CCL magazine. These incidents have 
aided greatly in the development of realistic scenarios to predict ship survivability. 

C. PREVIOUS WORK 

Bowman and Lee [Ref. 1] developed a one dimensional heat transfer model to 
predict the time-temperature profiles in the missile canisters of the MK 41 VLS system for 
Df^Cr-51 FLT I, II and IIA ships. Additionally, this model was used to evaluate the 
effectiveness of high temperature thermal insulation to reduce the heat transfer rates. 

Their scenarios used a fire in an adjacent compartment, which increased temperatures 
from ambient conditions to 2000°F in a period of five minutes, and then maintained that 
temperature for the duration of the fire. 

Mansfield [Ref 3] conducted an analysis of ordnance heating in steel walled 
compartments. He developed a computer algorithm that introduces heat into the fire 
compartment according to any selected fire size, and follows the heat transfer into the 
magazine and other compartments that surround the fire compartment. The output of this 
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algorithm is the calculation of the temperature of the fluid medium in the compartment and 
the surface temperature of the surrounding bulkheads, deck and overhead. It does not, 
however, take into account any varying geometry within the magazine. 

D. OBJECTIVES 

The purpose of this study is to model a particular CCL within the module, 
specifically the CCL located in the first row along the centerline and corner, for a DDG-51 
FLT-IIA ship. This one dimensional lumped parameter model describes the time- 
temperature profile in a Tomahawk missile to predict the time it takes to reach a critical 
value or “cook-ofF’ temperature of propellant at the geometric center of the missile. 
Different fire intensities, temperatures and heat fluxes will be used to provide a wide 
baseline for comparison to real world scenarios. 
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11. MODEL 


A. FORMULATION 

The computer algorithm modeling the transfer of heat from a compartment on fire 
to a missile contained in a Concentric Canister Launcher (CCL) magazine is based on the 
electrical circuit analogy for heat transfer. Just as an electrical circuit is associated with 
the conduction of electricity, a thermal circuit may be associated with the transfer of heat 
[Ref 4], Therefore, the first step in the creation of the model is to develop a thermal 
network describing the heat transfer from the space to the missile taking into account the 
geometry of the CCL magazine. Time dependence is accounted for in the network by 
assigning a capacitance for the solid material of the bulkhead, canister walls, and missile. 
The input parameter to the algorithm is the rate of temperature increase on the outer 
surface of the bulkhead from ambient conditions to some specified steady value based on 
the fire scenario analyzed. Thermal conduction, convection and surface radiation are all 
accounted for in the model. Each one will be discussed along with the associated 
correlations, thermophysical properties, and geometry associated with its position in the 
thermal resistance network. The computer code for the heat transfer model is written in 
MATLAB 4.2.C. A schematic of the thermal resistance network is given in Figure 4. 

1. Conduction 

Conduction takes place in four locations in the resistance network. The first is the 
bulkhead which is modeled as a plane wall constructed of mild steel, AISI 1020. The 
thermophysical property data for this material were obtained from Touloukian’s, 
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“Thermophysical Properties of Matter” [Ref 6], The data were curve fitted using 
MATLAB to form a polynomial to develop a smooth curve that “best fits” the property 
data and can be evaluated easily by the algorithm. The properties required for the 
algorithm are thermal conductivity (kj,,), specific heat (Cj,), and density (Ps,). Each were 
evaluated at the average temperature between the two surfaces of the bulkhead. The 
thermal resistance for conduction through the bulkhead is given by the equation 


^^cond 


Ax 


( 1 ) 


where (Ax) is the bulkhead thickness and (A,) is the area normal to heat transfer. 
Therefore, the thermal energy transferred through the solid due to conduction is described 
by the relationship; 


^cond 


T -T 

i 1 jj 


R1 


cond 


( 2 ) 


The capacitance assigned to the bulkhead to account for time dependence is expressed in 
the form of the rate of increase in stored energy, 


^stored 


‘ AT 

av 

“aT 


( 3 ) 
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Figure 4. Thermal Resistance Network of the CCL Magazine. 
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where the average temperature is used to describe the temperature profile through the 
solid. This model assumes temperature profiles to be linear at any instant through solid 
material. The capacitance equation for stored energy remains the same for the cylinder 
walls and missile, except in the thermophysical properties associated with the different 
material used in their construction and in the calculation of the volume, which takes into 
account their different geometric shapes. 

Conduction through the outer and inner cylinder walls is evaluated similar to the 
bulkhead, except that it is modeled in cylindrical coordinates and through walls 
constructed of titanium. Their thermophysical properties are again obtained from 
Touloukian, and are curve fitted with MATLAB to form a polynomial for use by the 
algorithm [Ref. 6], The thermal resistance equations for conduction in the outer and inner 
cylinder wall are; 
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cond 
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The missile’s thermal resistance is constructed in the same manner as above except the 
missile thermophysical properties are assumed to be aluminum and uniform throughout. 

2. Convection 

Natural convection occurs in three areas in the resistance network. The first is the 
space between the bulkhead and the outer surface of the CCL. The second is the annulus 
between the inner and outer cylinders, which function as the gas management system 
during missile launch. The third is the space between the inner cylinder and the surface of 
the missile. The fluid medium of each area is assumed to be air and is evaluated at a mean 
boundary layer temperature, termed the film temperature. The properties required for 
computing the convective resistance are thermal conductivity (k^ij), kinematic viscosity 
(v), volumetric thermal expansion coefficient (P), and Prandtl number (Pr). The thermal 
resistance for convection for each air space is given by the equation 


R = 


conv 


1 

h><A 


( 6 ) 


where the calculation of the average convection heat transfer coefficient ( /; ) is based 
upon the Nusselt number ( ), thermophysical properties of air, and either the height 

(L) of the CCL or the bulkhead. 

_ FxA 

h= " (7) 
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The thermophysical properties of air were curve fitted from tables in Reference 4 for use 
by the computer algorithm. The Nusselt number is a dimension less heat transfer 
coefficient and provides a measure of the convection heat transfer occurring at the surface. 
Empirical correlations for the Nusselt number have been developed based on the geometry 
of each enclosed air space. 

The first air space Nusselt numbers are calculated with respect to the bulkhead and 
the CCL’s free convection boundary layer. The bulkhead’s Nusselt number is determined 
by using a correlation for a vertical plate developed by Churchill and Chu that may be 
applied over an entire range of Rayleigh numbers [Ref 4]. 


( 



\ 


0.825 


0.387/?a/ 


1J 0.492] 

9 ' 
16 

[ [ Pr j 



( 8 ) 


The same equation may be applied to the outer surface of the CCL if the boundary layer 
thickness is much less than the diameter of the cylinder. To check for this condition the 
following criterion must be met [Ref 4]. 


d 



can 


> 



(9) 
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If this inequality is not met, then curvature effects are appreciable and must be accounted 
for using the following correlation by Le Fevre and Ede [Ref 5]; 



IGr^Pr'^ U^ 4(272 +3 5Pr)L 
5(20 +2\Pr) J 3 5(64 +62Pr)d 


( 10 ) 


The Grashof (Gr) and Rayleigh (Ra) numbers used throughout these correlations are given 
by: 




( 11 ) 


Ra.^Grfr 


( 12 ) 


After each Nusselt correlation is computed, the average convective heat transfer 
coefficient is calculated from Equation 7. Assuming both surfaces act as external free 
convection flows, resistances are calculated with respect to the boundary layers along the 
bulkhead and CCL’s surface using Equation 6. These resistances are added in series to 
form the total convective resistance between the bulkhead and the outer surface of the 
CCL. A schematic of the resistances is given in Figure 5. 
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Figure 5. 1st Air Space Resistance Network. 


The final two air spaces are circular annuli, and the Nusselt correlation used was 
developed by Churchill and Ozoe [Ref. 7]. 
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In the two circular annuli only one resistance is calculated to describe the 
resistance across the enclosed space. The resistance and average convective heat transfer 
coefficient is again calculated using Equation 6 and 7, respectively. A schematic of the 
resistances across the annuli is given in Figure 6. 



Figure 6. Resistance Network Across Annuli. 


3. Radiation 

Radiation affects heat transfer in the same air spaces as convection. Therefore, the 
thermal resistance for radiation and convection act in parallel, and can be combined to 
obtain a single effective surface resistance [Ref 4], 
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( 15 ) 
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conv 
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rad 


The total resistance to radiation exchange between the surfaces is comprised of the two 
surface resistances and the geometrical resistance. 


Rr 


tot 


-^ 23^2 ^ 3^3 


( 16 ) 


The radiative property used in the resistance calculations is emissivity (e). The emissivity 
of the different materials is obtained from Touloukian [Ref 6], and curve fitted to form a 
polynomial for use in the computer algorithm. The view factor (Fjj) is explained later. 
Therefore, the net radiation exchange between surfaces is expressed as; 


^o=- 


oiT'-T') 


Rr 


tot 


( 17 ) 


However, this equation is non-linear to the fourth power, and the overall radiative 
resistance (R^aa) must be solved, so the heat rate due to radiation takes the form; 


^rad 


^rad 


( 18 ) 
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and can be used in conjunction with the convective resistance. This is accomplished using 
two different methods. The first is to linearize the equation by expanding in a Taylor 
series, and the second is to solve it explicitly. Using a Taylor series expansion on 
Equation 17 takes the form; 




dT)r. 


(r-7;)+0[(r-r,)"] 


( 19 ) 


and after linearizing, the resistance due to radiation becomes 


<lo ^ 4or3 


Kad Rr 


tot 


( 20 ) 


To solve for the overall radiative resistance explicitly. Equations 17 and 18 are set equal to 
each other and solved algebraically. The overall resistance due to radiation is then 
expressed as: 


1 

^rad 





( 21 ) 


Note that the resistance is temperature dependent. These two processes are used for all 
radiative resistances in the thermal resistance network, and will be compared in the final 
results. The difference in the geometric shapes of the enclosures are taken into account by 
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calculating the fraction of the radiation leaving one surface and intercepted by another. 
This is termed the view factor (Fjj). In the first air space, the view factor is calculated as 
an infinite plane to a row of cylinders [Ref 9], 


I / 

F^^=\-{\-D^)^ +D arctan 

V 
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( 22 ) 


where 



(23) 


In the final two air spaces, the view factors are calculated as the interior of an outer right 
circular cylinder of finite length to the exterior of a inner right circular coaxial cylinder of 
finite length [Ref 10], 
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where 
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( 25 ) 


H= 



R= 


( 26 ) 


The thermal energy transferred through the air space is described by the relationship 


< 1 = 




( 27 ) 


where (R,) is the total resistance due to radiation and convection, Equation 15. 

The radiative resistance network thus far models a CCL on the first row near the 
centerline of the ship. As mentioned in the objectives, the CCL in the comer of the 
magazine is also modeled, accounting for the radiation effects of the adjoining bulkhead. 
The adjoining bulkhead is modeled as a reradiating surface. Equation 16 is modified to 
account for the radiation leaving the effected bulkhead, reflected off the adjoining 
bulkhead and absorbed by the surface of the CCL, and takes the form; 
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The view factor from the reradiating bulkhead to the canister (F^j) is the same as 
the view factor from the effected bulkhead to the canister (F 23 ), Equation 22. However, a 
new view factor (Fjr) to determine the amount of radiation received by the adjoining 
bulkhead is calculated as perpendicular plates with a common edge [Ref 4], 
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(29) 


These equations are solved as previously discussed. A schematic of the surface radiation 


resistances is given in Figure 7. 



Figure 7. 1 st Air Space Resistance Network (Reradiating Surface). 
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4. Finite Difference Equations 

The capacitance and transferred thermal energy from conduction and the effects of 
convection and radiation are combined to form a set of explicit finite difference equations 
to solve for the temperature at each node in the thermal resistance network. Each is 
solved simultaneously at each instant in time, which is defined by the time step specified in 
the algorithm. These finite difference equations were derived by conducting an energy 
balance about each temperature node. Because this is a time dependent model due to the 
capacitance in the solid material of the bulkhead, canister walls and missile, the thermal 
energy stored in the solid material is taken into account in the energy balance. 

K^^stored ( 30 ) 


All heat flow is assumed to flow into the node, and that each node in the resistance 
network is initially at ambient conditions (300 K). The result of the energy balance yields 
equations for the nodal temperature at a time increment. 


7’,(p+l)=2 


Ar 

pcAAx 


T,(P)-T^(P) 


At 

pc4Ax 


T,(p)-TJp)\ f T,(pyT,(p) 


- T,(p+l) 



21 



A/ 

' T,(p)-T,(p)' 

At 


+ 


pcAAx 

^ ^^cond J 

pcALx 

R2, 

1 


2 J 




( 33 ) 


At 

' W-T,(p) 

At 

' Te(p)-T,(p)\ ^ 


pcAAx 

R2, 

\ ^ ) 

pcAAx 

^^cond j 

2 J 




( 34 ) 


At 

' T,(p)-T,(p)' 

At 

+ 
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J 
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At 

' W-T^(p)^ 

At 



' T^^p)^T,(p)] 

pcAAx 

J 

pcAAx 

^‘^cond ) 


2 j 


7’8(P+1) 


( 36 ) 


7’,(p+l)=2 


At 

' w- 

+ 

T^(p)^T,(p)]' 

pcAAx 

^ ^^cond J 


2 J 




( 37 ) 


B. LIMITATIONS 

In the construction of this model, a high level of accuracy was not anticipated due 
to its one dimensionality. However, the goal is to develop a model that achieves 
moderately accurate solutions to provide baseline data on the cook-off of Tomahawk 
missiles contained in a CCL magazine. The primary limitations and associated errors with 
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this model are in the finite-differencing method, thermophysical property data, geometric 
assumptions, and simplification of the heat transfer process. 

1. Finite Difference Method 

The finite difference solution used is the explicit method. In this solution method, 
the unknown nodal temperatures for the new (p+1) time are determined from known 
nodal temperatures at the previous (p) time. The accuracy of this method depends on the 
distance between nodes (Ax) and the time step (At), as these values decrease the accuracy 
of the solution increases. The distance between nodes is chosen based on the geometry of 
the model, and the time step is chosen based on the stability requirements. The diflSculty 
with the explicit method is that it is not unconditionally stable. This condition is 
characterized by oscillations which diverge from a steady value instead of converging as 
anticipated in transient conduction problems. To determine the time step at which the 
algorithm used in the model is stable, the following criterion for a one dimensional node is 
solved for At. 


ctAr _ 1 
(Ax)2 2 


( 38 ) 


The limiting time step in this model is governed by the inner cylinder wall based on its 
values of wall thickness (Ax) and thermal difiiisivity (a). Therefore, maximum allowable 
time step is 0.37 seconds. The more the time step is reduced below 0.37 seconds, the 
more accurate the solution. However, the number of iterations to achieve the desired time 
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span increases rapidly and makes it computationally intensive. Therefore, 0.37 seconds is 
used in all calculations to optimize the computing time. 

2. Property Data 

The thermophysical properties obtained from Touloukian [Ref. 6], which are 
entering arguments throughout the algorithm, are a large source of error. Not only is 
there a error associated with curve fitting of the data, but also in the experimental 
procedure in obtaining the properties. The error reported by Touloukian [Ref 6] in the 
experimental procedures used to obtain the thermophysical property data of the substances 
used in the CCL are listed in Table 1. 


EXPERIMENTAL ERROR 


STEEL 

TITANIUM 

AIR 

Specific Heat (cp) 

4% 

0.2% 

n/a 

Thermal Conductivity (k) 

1% 

5-15% 

1-5% 

Emissivity (e) 

2% 

8% 

n/a 

Thermal Dififiisivity (a) 

5-10% 

10% 

n/a 

Kinematic Viscosity (v) 

n/a 

n/a 

2% 


Table 1. Experimental Error in the Measurement of Thermophysical Properties. 


In order for the algorithm to evaluate the properties of the metals and air, a curve 
fit or regression was conducted on the data. In this procedure, a curve was developed that 
‘best fits” the data but does not necessarily pass through the data points. The “best fit” or 
least squares curve fitting seeks to minimize the sum of the squared error. The error is the 
distance between the polynomial curve and the actual data point. Squaring this distance at 


24 

























each data point and adding the squared distances together is the sum of the squared error. 
Therefore, the curve generated by the polynomial is the curve that makes the sum of the 
errors as small as it can be, “best fit” [Ref 8], The error in relation to the data and the 
polynomial generated by the curve fit was less than ten percent in all cases, which is within 
the range of the experimental data error. 

3. Geometry 

In the modeling of the CCL each geometric shape is considered to be solid 
polished material with no attachments. Obvious architecture with respect to a functioning 
weapon system have been neglected such as wire ways, electronic packaging, supporting 
fi’ame work, etc. 

Several assumptions have been made with respect to radiation. The metal on all 
surfaces throughout the model are assumed to be unpainted and polished. Thus, the 
emissivity is assumed to be much higher than that of a painted surface. Additionally, the 
model accounts only for the heat emitted from the bulkhead to the modeled CCL and 
neglects radiation effects fi-om the surrounding canisters. Finally, radiation with respect to 
the other walls in the room on fire to the CCL bulkhead is neglected, as well as, the 
convective resistance on the outer surface of the bulkhead. 

The final assumption in the model is that the missile has the properties of aluminum 
and is solid throughout. The missile electronics and cabling are neglected. Aluminum was 
chosen because it is used consistently throughout the missile. Its property data would 
provide a good assumption as to how the missile reacts to the heat energy absorbed and 
provide a temperature profile at the center of the missile. 
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III. FIRE SCENARIOS 


To develop accurate time-temperature profiles describing the heat-up rate of the 
CCL, credible fire scenarios are identified based on the type of combustible material 
fueling the fire and its probable intensity. Four scenarios are analyzed by the CCL heat 
transfer model to achieve a good understanding of how the CCL’s heat-up rate will differ 
when exposed to different types of fires. 

The first scenario analyzed consisted of a compartment fire involving common 
shipboard fuels (JP-5, JP-4, F-76, etc). The time-temperature profile is characterized by a 
increase in temperature from ambient conditions to 2000"F (1366 K) over a period of five 
minutes and is then maintained at the maximum temperature for the duration of the fire 
[Ref 1]. This scenario is typically used to evaluate the structural integrity of the effected 
space both during and after the fire [Ref 11]. 

The second fire scenario is typically used in modeling fires fiieled by common 
combustible material, such as a berthing compartment fire. The time-temperature profile 
is characterized by a increase in temperature from ambient conditions to 1700"F (1200 K) 
over a period of fifteen minutes, and then maintains a steady temperature rise of 84"F/hour 
(46 K/hour) [Ref 11]. 

A fire fueled by a missile’s residual solid propellant, as experienced by the USS 
STARK (FFG 31), is the most intense and is the next scenario analyzed. Propellant fires 
are normally short lived but can reach temperatures in excess of 3000T (1922 K). For the 
purpose of this scenario, it is assumed that maximum temperature is reached in thirty 
seconds, and that the temperature is then maintained until the self-ignition temperature of 
the missile’s solid propellant is reached. 
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The last scenario analyzed was developed by the U. S. Navy’s Insensitive 
Munitions (IM) Program to characterize the thermal response of munitions exposed to 
slow cook-off conditions. The slow cook-off condition is described by a temperature rise 
of6T/hour(3.33 K/hour) [Ref 1], 

The IM Program also developed a fire scenario to describe fast cook-off, that 
immerses the missile in a JP-5 pool fire in which the average temperature must be at or 
above 1600T [Ref 1], This scenario is most likely on ordnance stored on a flight deck 
rather than a missile contained in a CCL magazine. Therefore, the fast cook-off scenario 
is not analyzed. 

The test data from the Insensitive Munitions program shows that the missile’s solid 
propellant motor (MK 106 Booster) and explosives will self-ignite between 300 and 400°F 
(422 and 478 K), and the liquid propellant motor, fueled by JP-10 or RJ-4, can self-ignite 
at temperatures as low as 460“F (511 K). Depending on the rate of temperature increase, 
self-ignition of the energetic material can range from a slow burn to an immediate 
deflagration. [Ref 1] 

Table 2 summarizes the aforementioned fire scenarios and assigns a scenario 
number. This number is used to refer to the various fire scenarios throughout the rest of 
the text. 
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FIRE SCENARIOS 


MAXIMUM 
TEMPERATURE (K) 

RATE OF TEMP 
INCREASE 

1366 

3.55 K/second 

n/a 

80K/niin< 1186K 

0.77 K/min > 1186 K 

1922 

54.07 K/second 

n/a 

3.33 K/hour 


Table 2. Fire Scenarios. 
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IV. RESULTS 


As mentioned in the objectives, this study models two CCLs within the ship’s 
weapon module. The first is located in the first row along the centerline of the ship, and 
the second is the corner CCL in which the adjacent wall is considered a reradiating 
surface. The four fire scenarios mentioned in the previous chapter are analyzed for each 
CCL location. Additionally, each scenario is conducted twice. The first solves for surface 
radiative resistance explicitly and the second with a Taylor series expansion. After the 
scenario number, the explicit method and Taylor series expansion methods will be denoted 
by the letter (A) and (B), respectively. 

The temperature range in which the missile’s components, specifically the warhead 
and solid propellant motor, will self-ignite is between 422 and 478 K. The liquid 
propellant will self-ignite as low as 511 K. The times to reach these temperatures for the 
center and comer CCL are listed in Table 3 and 4, respectively. Time-temperature profiles 
are plotted in Figures 8 through 15 for the center CCL and Figures 16 through 23 for the 
comer CCL. 

The two methods used to determine the surface radiation, described in Chapter II, 
are compared because of the large differences experienced in each scenario’s time- 
temperature profiles. The Taylor series expansion produces two types of error, tmncation 
and round-off. The truncation error is due to the approximation formula. In general, 
increasing the order of the Taylor series reduces the tmncation error but increases the 
round-off error. For this model both types of errors are probably high for two reasons: 
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the expansion was only carried out to the first derivative to reduce the amount of 
computations, and the round off error is high due to the curve fitted thermophysical 
properties used throughout the model. Solving for the surface radiative resistance 
explicitly eliminates the truncation error, however, it is still affected by the round off error. 
Solving for the surface radiative resistances with a Taylor series expansion produces a 
more conservative time-temperature profile, and is the recommended method. 


32 




TIME TO COOK-OFF 

SCENARIO # 

TIME TO REACH 
511 K (min) 

TIME TO REACH 
422-478 K (min) 

ASSOCIATED 

FIGURE 

lA 

204.42 

141.06-181.02 

8 

IB 

152.26 

108.48-136.32 

9 

2A 

254.40 

187.52-230.96 

10 

2B 

198.68 

147.43-180.73 

11 

3A 

64.85 

49.71-59.28 

12 

3B 

53.25 

42.23-49.25 

13 

4A 

inf 

inf 

14 

4B 

inf 

inf 

15 


Table 3. Cook-Off Times for Scenarios (Center CCL). 


TIME TO COOK-OFF 

SCENARIO # 

TIME TO REACH 
511 K (min) 

TIME TO REACH 
422-478 K (min) 

ASSOCIATED 

FIGURE 

lA 

199.88 

137.76-176.95 

16 

IB 

149.59 

106.43-133.88 

17 

2A 

250.04 

183.92-226.88 

18 

2B 

195.83 

145.02-178.03 

19 

3A 

63.50 

48.67-58,06 

20 

3B 

52.42 

41.56-48.48 

21 

4A 

inf 

inf 

22 

4B 

inf 

inf 

23 


Table 4. Cook-Off Times for Scenarios (Comer CCL). 
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Figure 8. Scenario lA-Missile Heat-Up Rate. 
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Figure 9. Scenario IB-Missile Heat-Up Rate. 
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Figure 10. Scenario 2A-Missile Heat-Up Rate. 
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Figure 13. Scenario 3B-Missile Heat-Up Rate. 
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Figure 16. Scenario 1 A-Missile Heat-Up Rate w/ Reradiating Surface. 
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Figure 20. Scenario 3 A-Missile Heat-Up Rate w/ Reradiating Surface. 
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Figure 21. Scenario 3B-Missile Heat-Up Rate w/ Reradiating Surface. 
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Figure 23. Scenario 4B-Missile Heat-Up Rate w/ Reradiating Surface. 
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V. CONCLUSIONS 


The fire scenarios analyzed by the CCL heat transfer model in this study showed 
that the CCL can withstand an intense conflagration if the bulkheads surrounding the 
magazine remain intact. In the most intense fire situation, Scenario 3, ship’s personnel 
have approximately forty-five minutes to gain control of the fire or cool the weapons in 
the magazine prior to cook-off. Conceivably, automated fire protection equipment could 
increase the amount of time firefighters have to control the fire without the immediate 
hazard of cook-off 

In its current stage of development, the CCL’s installed fire protection features 
have not been designed. However, the uniqueness of the CCL module will require some 
changes, as compared to the fire protection equipment currently installed in the MK 41 
VLS. The MK 41 VLS fire protection system consist of a deluge system installed in the 
canister designed to cool the missile in the event of an inadvertent motor ignition or 
extreme over temperature situation, and a sprinkling system designed to provide cooling 
external to its canisters should it sense a rapid temperature rise [Ref 1]. Because the CCL 
acts both as a launcher and transport device, introducing a fire protection system into the 
CCL will hamper its primary advantage of being a totally integrated package. 

A possible alternative to the current deluge system would be a similar system of 
two or more mounted spray nozzles directed to impinge the canister at the top outer 
surface of the CCL providing a cooling film around the entire canister but not wetting the 
missile. In addition to the outer deluge system, a magazine sprinkler should also be 
installed to cool the entire magazine should temperatures increase too rapidly. The 
advantages of keeping the fire protection system outside of the CCL are two fold; 
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• If there is an inadvertent sprinkling of the magazine or activation of the 
deluge system, the missile is not effected by the water and the costly 
process of sending the weapon to a rework facility for repairs and 
recertification or decertification and dismantling is avoided. It is possible 
that after a incident of this type onboard data transmission, continuity and 
accuracy checks can be conducted to ve^rify the missile is fully functional 
after the module and electrical connections are cleaned. 

• The second advantage, and most important, is to maintain the ease of the 
installation and removal process from the ship’s weapon module. Installing 
piping connections for a deluge system, which of course will have to meet 
high shock criteria, will make a relatively easy installation/removal 
operation much more cumbersome. 

The disadvantage of water systems is the generation of steam in intense fires 
resulting in a steam explosion. There are other types of systems, such as inert gas fire 
protection systems, which may be worthy of investigation. But when dealing with intense 
fires with long durations, as in the USS STARK, the availability and cooling effects of 
seawater proved to be effective. Investigating the installation of a relief system should 
steam accumulate in the module during a fire is recommended. 

Another characteristic of missile solid propellent, brought to our attention by 
NAVAIR, China Lake, which warrants mentioning, is that prior to reaching self-ignition 
temperature, the propellant begins to off-gas some of its volatile materials [Ref 11]. The 
temperature at which this occurs is yet to be identified. From the standpoint of any 
exposure to an over-temperature situation, a missile may be damaged even if cook-off did 
not occur. However, an interesting side effect was noticed that if the solid propellant was 
remixed after the off gassing, it burned longer and with less thermal energy. Follow on 
research is recommended to determine at what temperature the solid propellent begins to 
off-gas. 
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In general, this study will provide CCL project managers and designers with 
baseline time-temperature profiles to make decisions concerning the construction and 
survivability of the CCL launching system, especially with respect to installed fire 
protection systems. Additionally, it will provide ship’s personnel with a conservative 
missile cook-ofF time based on the type of fire, so that the damage control situation can be 
assessed and priorities set. 
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VI. RECOMMENDATIONS 


In any continuation of this study the following are recommended: 

• The CCLs modeled are sized for Tomahawk missiles and are the same 
throughout the ship’s weapon module. Future studies should model the 
ship’s weapon module for an actual deployment weapon configuration in 
which different types of CCLs carry a variety of weapons in support of all 
the ship’s mission areas. 

• Build a user interface at the MATLAB command window to easily vary the 
fire scenario. Currently, to vary the fire scenarios the computer code must 
be modified. 

• Account for surface radiation from other CCLs in the magazine. 

• Utilize temperature data from fire scenarios currently being conducted in 
the ex-US S SHAD WELL. This should provide more accurate data and a 
larger variety of fire curves to be evaluated by the CCL Fire Model. 
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APPENDIX A. CCL-FIRE PROGRAM CODE (2A-CENTER) 

This MATLAB code solves for radiation explicitly as described in the radiation 


section of Chapter II. This code is set up for the center CCL and scenario 2. 


format bank 
clear 

% DIMENSIONS FOR CANISTERS & SURROUNDING VLS WALL. 


w=l; 

% 

Lw=l; 

% 

Lc=l; 

% 

Lm=l; 

% 

delx=0.0048; 

% 

r=0.00159; 

% 

ro=0.2586; 

% 

rl=0.2635; 

% 

r2=0.2667; 

% 

r3-0.3492; 

% 

r4=0.3571; 

% 

s=0.8382; 

% 

g=9.81; 

% 

sig=5.67e-008; 

% 

Al=w*Lw; 

% 

A2=A1; 


A3=2*pi*r4*Lc; 

% 

A4=2*pi*r3*Lc; 

% 

A5=2*pi*r2*Lc; 

% 

A6=2*pi*rl*Lc; 

% 

A7=2*pi*ro*Lc; 

% 

Vl=deix*w*Lw; 


V2=pi*Lc*(r4^2-r3^2); 

V3=pi*Lc*(r2'^2-rl^2); 

V4=pi*ro'^2*Lm; 



[m] Width of for and aft bulkhead. 

[m] Height of wall surrounding CCL module. 

[m] Height of inner/outer canister wall. 

[m] Length of missile. 

[m] Thickness of wall surrounding CCL module, 
[m] Near center of missile (1/16 in). 

[m] Radius of missile. 

[m] Inside radius of inner cylinder wall. 

[m] Outside radius of inner cylinder wall. 

[m] Inside radius of outer cylinder wall. 

[m] Outside radius of outer cylinder wall. 

[m] Distance between cylinder CL. 

[m/s'^2] Gravitational constant. 

[W/m^2*KM] Stefan-Boltzmann constant. 
[m^2]Area of VLS wall perp to dir of heat flux. 

[m^2] Outer surface area of 1st cylinder. 

[m'^2] Inner surface area of 1st cylinder. 

[m'^2] Outer surface area of 2nd cylinder. 

[m^2] Inner surface area of 2nd cylinder. 

[m'^2] Surface area of missile. 


% [K] INITIAL AMBIENT CONDITIONS TO START ITERATION. 

T2(l)=300; 

T3(l)=300; 
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T4(l)=300; 

T5(l)=300; 

T6(l)=300; 

T7(l)=300; 

T8(l)=300; 


% DEFINES TIME STEP AND NUMBER OF ITERATIONS 

S=50000; % Number of iterations. 

delt=0.36; % [s] Time step. 


% BEGINS LOOP TO SOLVE FOR ALL TEMPS IN RESISTANCE NETWORK. 


for p=l;S 
pl=p 


% [s] Time. 
t(p)=(delt*(p-l)); 

dT(p)=t(p) * 1.33 3 3; 
dTA(p)=t(p)*0.0128; 
Tl(p)=300+dT(p); 
TlA(p)-1200+dTA(p); 

ifTl(p)<=1200 

Tl(p)=Tl(p); 

else 

Tl(p)=TlA(p); 

end 

t(p+lHdelt*(p)); 
dT(p+l)=t(p+l)* 1.3333; 
dTA(p+l)=t(p+l)*0.0128; 
Tl(p+l)=300+dT(p+l); 
TlA(p+l)=1200+dTA(p+l); 

ifTl(p+l)<=1200 

Tl(p+l)=Tl(p+l); 

else 

Tl(p+l)=TlA(p+l); 

end 
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%%% SOLVES SECTION 1 [T2] %%% 


% [K] Film temperature to evaluate properties of mild steel 
Tfl(pHTl(p)+T2(p)y2; 

% [W/m*K] Thermal conductivity of bulkhead. 
kwl(p)=abs((1.8869e-005*Tfl(p)^2)-(6.9896e-002*Tfl(p))+83.438); 

% [J/kg*K] Specific Heat of bulkhead. 

cp 1 (p)=abs(((-2.2828e-006)*Tfl(p)^3)+((4.6812e-003)*Tfl (p)^2)-(2.3469*Tfl (p))+800. 
99); 

% [m^2/s] Thermal diffusivity of bulkhead. 

alfl(p)=abs(((4.0118e-014)*Tfl(p)^3)-((7.9337e-011)*Tfl(p)^2)+(3.3733e-008*Tfl(p)) 

+1.0388e-005); 

% [kg/m^3] Density of bulkhead, 
rho 1 (p)=kwl (p)/(cp 1 (p)*alfl (p)); 

% Solves for the conductive resistance in bulkhead. 

Rcl (p)=delx/(kwl (p)* A2); 

% [K] Film temperature to evaluate properties in 1st air space. 

Tf2(p)=(T2(p)+T3(p))/2; 

% [W/m*K] Thermal conductivity of air. 

ka2(p)=abs((4.9198e-011*T£2(p)^3)-(1.608e-007*Tf2(p)^2)+(2.014e-004*T£2(p))-0.020 

796); 

% [m^2/s] Kinematic viscosity of air. 

nu2(p)=abs(((7.786e-005*Tf2(p)^2)+(4.4053e-002*Tf2(p))-2.4972)* 1 .Oe-006); 

% Prandtl #. 

Pr2(p)=abs((-2.5065e-011 *T£2(p)^3)+(8.4944e-008*T£2(p)^2)-(l .0417e-004*Tf2(p))+0. 
74553); 


% [K^-1] Exp. Coefficient of air. 

B2(p)=l/Tf2(p); 

% Emissivity of mild steel. 

E2(p)=abs((-9.8416e-008*T2(p)^2)+(2.6912e-004*T2(p))+0.73699); 
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% Emissivity of titanium. 

E3(p)=abs((-5.0412e-008*T3(p)^2)+(2.4949e-004*T3(p))+0.051504); 

% Solves for Rayleigh and Girashof #'s WRT bulkhead inner surface. 
Gr2(pH(g*B2(p)*(T2(p)-Te(p))*(Lw-3))/(nu2(pr2)); 

Ra2^)=Gr2(p)*Pr2(p); 

% Correlation to solve for the average Nusselt # over a 
% vertical flat plat. Recommended by Churchill and Chu, it is 
% good over entire range of Rayleigh #'s. 

Nu2(p)-(0.825+((0.387*(Ra2(p)^(l/6)))/((l+((0.492/Pr2(p))^(9/16)))^(8/27))))^2; 


% Solves for convective resistance using the average convection 
% heat transfer coefficient over the bulkhead inner surface. 
h2(p)=Nu2(p)*ka2(p)/Lw; 

Rco 1 a(p)= 1 /(h2(p) * A2); 


% Solves for Rayleigh and Grashof #'s WRT cylinder outer surface. 
Gr3(pHg*B2(p)*(Tf2(p)-T3(p))*(Lc-3))/(nu2(pr2); 

Ra3(p)=Gr3(p)*Pr2(p); 

% IF-THEN LOOP TO DETERMINE IF THE OUTER WALL OF THE CYLINDER 
% CAN BE APPROXIMATED AS A PLANE WALL. IF A PLANE WALL 
% CHURCHILL 

% AND CHU'S CORRELATION FOR A VERTICAL SURFACE IS USED. IF 
% THERE ARE APPRECIABLE CURVATURE EFFECTS A CORRELATION 
% DEVELOPED 

% BY LE FEVRE AND EDE BASED ON CYLINDER HEIGHT. 


if Gr3(p)==0 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)^2)/(5*(20+21*Pr2(p))))^0.25)+((4*Lc*(272+315*Pr2 

(p)))/(35*2*r4*(64+63*Pr2(p)))); 

else 

if ((2*r4)/Lc)>=3 5/(Gr3(p)^l/4) 

Nu3(pH0.825+((0.387*(Ra3(p)^(l/6)))/((l+((0.492/Pr2(p))^(9/16)))^(8/27))))^2; 

else 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)^2)/(5*(20+21*Pr2(p))))^.25)+((4*Lc*(272+315*Pr2 

(p)))/(35*2*r4*(64+63*Pr2(p)))); 

end 

end 
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% Solves for convective resistance using the average convection 
% heat transfer coefficient over outer surface of the cylinder. 
h3 (p)=(Nu3 (p)*ka2(p))/Lc; 

Rco lb(p)= 1 /(h3 (p)* A3); 

% Adds convective resistances in series. 

Rco 1 (p)=Rco 1 a(p)+Rco 1 b(p); 

% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, infinite plane to a row of pipes [C-6]. 

D=(2*r4)/s; 

F23=l -((1 -D^2)^0.5)+(D*atan((( 1 -D^2)/D^2)^0.5)); 

% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

Rrl(p)=((l-E2(p))/(E2(p)*A2))+(l/(A2*F23))+((l-E3(p))/(E3(p)*A3)); 

Rradl(p)=l/((sig/Rrl(p))*((T2(p)+T3(p))*(T2(p)'^2+T3(p)^2))); 

% Solves for total resistance, [T3], and begins next iteration. 

Rt(p)=(Rco 1 (p)*Rrad 1 (p))/(Rrad 1 (p)+Rco 1 (p)); 

% Finite difference equation—explicit method. 

H1 (p)=delt/(rho 1 (p)*cp 1 (p)* V1); 

T2(p+l)=2*(((Hl(p)/Rcl(p))*(Tl(p)-T2(p)))+((Hl(p)/Rt(p))*(T3(p)-T2(p)))+((T2(p)+T 

l(p))/2))-Tl(p+l); 

ifT2(p+l)<300 

T2(p+1)=300; 

else 

T2(p+l)=T2(p+l); 

end 


%%% SOLVES SECTION 2 [T3] %%% 

% [K] Film temperature to evaluate properties of titanium. 

TD(p)KT3(p)+T4(p))/2; 

% [W/m*K] Thermal conductivity of canister wall. 

kc3(p)=abs((l .701 le-008*TD(p)^3)-(4.3316e-005*Tf3(p)^2)+(2.7279e-002*Tfi(p))+20 
358); 
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% [J/kg*K] Specific Heat of canister wall. 

cp3(p)=abs(((4.6480e-007)*Tf?(p)^3)-((1.4813e-003)*TO(p)^2)+(1.5387*T£3(p))+139.9 

3); 

% [m'^2/s] Thermal dififiisivity of canister wall. 

alf3(p)=abs(((-4.532e-015)*Tf?(p)^3)+((1.7281e-011)*Tf3(p)^2H1.929e-008*Tf3(p))+l 

.3594e-005); 

% [kg/m^3] Density of canister wall. 
rho3(p)=kc3(p)/(alf3(p)*cp3(p)); 

% Solves for the conductive resistance through the cylinder wall. 
Rc2(p)=log(r4/r3)/(2*pi*kc3(p)*Lc); 

% Finite difference equation-explicit method. 

H2(p)=delt/(rho3 (p) * cp3 (p)* V2); 

T3(p+l)=2*(((H2(p)/Rt(p))’^(T2(p)-T3(p)))+((H2(p)/Rc2(p))*(T4(p)-T3(p)))+((T3(p)+T 

4(p))/2))-T4(p); 


%%% SOLVES SECTION 3 [T4] %%% 

% [K] Film temperature to evaluate air properties in 1st annulus. 

Tf4(p)=(T4(p)+T5(p))/2; 

% [W/m*K] Thermal conductivity of air. 

ka4(p)=abs((4.9198e-011*Tf4(p)^3)-(1.608e-007*Tf4(p)^2)+(2.014e-004*Tf4(p))-0.020 

796); 

% [m'^2/s] Kinematic viscosity of air. 

nu4(p)=abs(((7.786e-005*Tf4(p)^2)+(4.4053e-002*Tf4(p))-2.4972)*1.0e-006); 

% Prandtl #. 

Pr4(p)=abs((-2.5065e-011 *Tf4(p)^3)+(8.4944e-008*Tf4(p)^2)-(l .0417e-004*Tf4(p))+0. 
74553); 

% [K^-1] Exp. Coefficient of air. 

B4(p)=l/Tf4(p); 


% Emissivity of titanium. 

E4(p)=abs((-5.0412e-008*T4(p)^2)+(2.4949e-004*T4(p))+0.051504); 

E5(p)=abs((-5.0412e-008*T5(p)'^2)+(2.4949e-004*T5(p))+0.051504); 
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% Solves for Rayleigh and Grashof #'s. 
Gr4(p)K(g*B4(p)*(T4(p)-T5(p))*(Lc^3))/(nu4(pr2)); 

Ra4(p)=Gr4(p)*Pr4(p); 

ifRa4(p)=0 

Ra4(p)=0.001; 

else 

Ra4(p)=Ra4(p); 

end 

% Nusselt correlation for circular annuli heated and cooled on the 
% vertical curved surface. Developed by de Vahl Davis and Thomas. 
fpr4(p)=(l+((0.5/Pr4(p))'^(9/16)))^(-16/9); 
Nu4(p)=0.364*((Ra4(p)*fpr4(p))^0.25)*((r3/r2)^0.5); 

% Solves for convective resistance using the average convection 
% heat transfer coefficient within the annuli. 
h4(p)=(Nu4(p)*ka4(p))/Lc; 

Rco2(p)=(A4+A5)/(h4(p)*A4*A5); 

% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, interior of outer right circular cylinder 
% of finite length to exterior of inner right circular coaxial 
% cylinder [C-87]. 

R=r3/r2; 

HH=Lc/r2; 

F45=(l/R)*(l-((HH'^2+R^2-l)/(4*HH))-l/pi*(acos((HH^2-R^2+l)/(HH^2+R^2-l))-(((sq 

rt(((HH^2+R^2+l)^2)-(2*(R^2))))/(2*HH))*acos((HH^2-R^2+l)/(R*(HH^2+R^2-l))))-( 

((HH^2-R^2+l)/(2*HH))*asin(l/R)))); 

% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

Rr2(p)=((l-E4(p))/(E4(p)*A4))+(l/(A4*F45))+((l-E5(p))/(E5(p)*A5)); 

Rrad2(p)=l/((sig/Rr2(p))*((T4(p)+T5(p))*(T4(p)^2+T5(p)^2))); 

% Solves for total resistance. 

Rt2(p)=((Rco2(p)*Rrad2(p))/(Rco2(p)+Rrad2(p))); 

% Finite difference equation-explicit method. 

T4(p+l)=2*(((H2(p)/Rc2(p))*(T3(p)-T4(p)))+((H2(p)/Rt2(p))*(T5(p)-T4(p)))+((T3(p)+ 

T4(p))/2))-T3(p+l); 
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ifT4(p+l)<300 

T4(p+1)=300; 

else 

T4(p+l)=T4(p+l); 

end 


%%% SOLVES SECTION 4 [T5] %%% 

% [K] Film temperature to evaluate properties of titanium. 

Tf5(p)=(T5(p)+T6(p))/2; 

% [W/m*K] Thermal conductivity of canister wall. 

kc5(p)=abs((1.7011e-008*T5(p)^3)-(4.3316e-005*Tf5(p)^2)+(2.7279e-002*Tf5(p))+20. 

358); 

% [J/kg*K] Specific Heat of canister wall. 

cp5(p)=abs(((4.6480e-007)*Tf5(p)'^3)-((1.4813e-003)*Tf5(p)^2)+(1.5387*Tf5(p))+139.9 

3); 

% [m^2/s] Thermal difilisivity of canister wall. 

alf5(p)=abs(((-4.532e-015)*Tf5(p)^3)+((l .7281e-01 l)*Tf5(p)^2)-(l .929e-008*Tf5(p))+l 
.3594e-005); 

% [kg/m^3] Density of canister wall. 
rho5(p)=kc5(p)/(alf5(p)*cp5(p)); 

% Solves for the conductive resistance through the cylinder wall. 
Rc3(p)=log(r2/rl)/(2*pi*kc5(p)*Lc); 

% Finite difference equation-explicit method. 

H3(p)=delt/(rho5(p)*cp5(p)*V3); 

T5(p+l)=2»(((H3(pyRt2(p))*(T4(p)-T5(p)))+((H3(p)/Rc3(p))*(T6(p)-T5(p)))+((T5(p)+ 

T6(p)V2))-T6(p); 


%%% SOLVES SECTION 5 [T6] %%% 

% [K] Film temperature to evaluate air properties in 2nd aimulus. 
Tf8(pHT6(p)+T7(p)y2; 


64 



% [W/m*K] Thermal conductivity. 

ka6(p)=abs((4.9198e-011*Tf6(p)^3)-(1.608e-007*Tf6(p)^2)+(2.014e-004*Tf6(p))-0.020 

796); 

% [m^2/s] Kinematic viscosity. 

nu6(p)-abs(((7.786e-005*Tf6(p)^2)+(4.4053e-002*Tf6(p))-2.4972)*1.0e-006); 

% Prandtl #. 

Pr6(p)=abs((-2.5065e-011 *Tf6(p)^3)+(8.4944e-008*Tf6(p)^2)-(l .0417e-004*Tf6(p))+0. 
74553); 

% [K'^-l] Exp. Coefficient of air. 

B6(p)=l/Tf6(p); 

% Emissivity of titanium. 

E6(p)=abs((-5.0412e-008*T6(p)^2)+(2.4949e-004*T6(p))+0.051504); 

E7(p)=abs((-5.0412e-008*T7(p)^2)+(2.4949e-004*T7(p))+0.051504); 

% Solves for Rayleigh and Grashof #'s. 
Gr5(p)=(g*B6(p)*(T6(p)-T7(p))*(Lc-3))/(nu6(p)-2); 

Ra5(p)=Gr5(p)*Pr6(p); 

ifRa5(p)=0 

Ra5(p)=0.001; 

else 

Ra5(p)=Ra5(p); 

end 

% Nusselt correlation for circular annuli heated and cooled on the 
% vertical curved surface. Developed by de Vahl Davis and Thomas. 
fpr5(p)=(l+((0.5/Pr6(p))^(9/16)))^(-16/9); 
Nu5(p)=0.364*((Ra5(p)*fpr5(p))^0.25)*((rl/ro)^0.5); 

% Solves for convective resistance using the average convection 
% heat transfer coefficient within the annuli. 
h5(p)=Nu5(p)*ka6(p)/Lc; 

Rco3(p)=(A6+A7)/(h5(p)*A6*A7); 

% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, interior of outer right circular cylinder 
% of finite length to exterior of inner right circular coaxial 
% cylinder [C-87]. 

Rl=rl/ro; 
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HHl=Lc/ro; 

F67=(l/Rl)*(l-((HH1^2+Rr2-l)/(4*HHl))-l/pi*(acos((HH1^2-Rr2+l)/(HHr2+Rr 

2-l)H((sqrt(((HH1^2+R1^2+l)^2)-(2*(R1^2))))/(2*HHl))*acos((HH1^2-R1^2+l)/(Rl* 

(HHr2+Rl'^2-l))))-(((HH1^2-Rl'^2+l)/(2*HHl))*asin(l/Rl)))); 

% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

Rr3(p)=((l-E6(p))/(E6(p)*A6))+(l/(A7*F67))+((l-E7(p))/(E7(p)*A7)); 

Rrad3(p)= l/((sig/Rr3(p))n(T6(p)+T7(p))*(T6(pr2+T7(pr2)^^ 

% Solves for total resistance. 

Rt3(p)=(Rco3(p)*Rrad3(p))/(Rco3(p)+Rrad3(p)); 

% Finite difference equation—explicit method. 

T6(p+l)=2*(((H3(p)/Rc3(p))*(T5(p)-T6(p)))+((ro(p)/Rt3(p))*(T7(p)-T6(p)))+((T 

T6(p))/2))-T5(p+l); 

ifT6(p+l)<300 

T6(p+1)=300; 

else 

T6(p+l)=T6(p+l); 

end 

% [K] Film temperature to evaluate properties of aluminum. 

Tf7(p)=(T7<p)+T8(p))/2; 

% [W/m*K] Thermal conductivity of missile. 
kc7(p)=237; 

% [J/kg*K] Specific Heat of missile. 
cp7(p)=903; 

% [m^2/s] Thermal diffusivity of missile. 
alf7(p)=97.1e-006; 

% [kg/m^3] Density of missile. 
rho7(p)-2702; 

% Solves for the conductive resistance through the cylinder wall. 
Rc4(p)=log(ro/r)/(2*pi*kc7(p)*Lc); 
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% Finite difference equation—explicit method. 

H4(p)=(delt/(rho7(p)*cp7(p)*V4)); 

T7(p+l)=2*(((H4(p)/Rt3(p))nT6(p)-T7(p)))+((H4(p)/Rc4(p))*(T8(p)-T7(p)))+((T7(^^^ 

T8(p))/2))-T8(p); 

T8(p+l)=2*(((H4(p)/Rc4(p))*(T7(p)-T8(p)))+((T7(p)+T8(p))/2))-T7(p+l); 

ifT8(p+l)<300 

T8(p+1)=300; 

else 

T8(p+l)=T8(p+l); 

end 

end 

t=t/60; 

plot(t,T8); grid; 
xlabel('Time (min)'); 
ylabel('Temperature (K)'); 
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APPENDIX B. CCL-FIRE PROGRAM CODE (2B-CENTER) 

This MATLAB code solve for radiation resistance using the Taylor series expansion as 


described in the radiation section of Chapter II. This code is set up for the center CCL and 

scenario 2. 

format bank 
clear 


% DIMENSIONS FOR CANISTERS & SURROUNDING VLS WALL. 


w=l; 

% 

Lw=l; 

% 

Lc=l; 

% 

Lm-1; 

% 

delx=0.0048; 

% 

r=0.00159; 

% 

ro=0.2586; 

% 

rl=0.2635; 

% 

r2=0.2667; 

% 

r3=0.3492; 

% 

r4=0.3571; 

% 

s=0.8382; 

% 

g=9.81; 

% 

sig=5.67e-008; 

% 

Al=w*Lw; 

% 

A2=A1; 


A3=2*pi*r4*Lc; 

% 

A4=2*pi*r3*Lc; 

% 

A5=2*pi*r2*Lc; 

% 

A6=2*pi*rl*Lc; 

% 

A7=2*pi*ro*Lc; 

% 

Vl=delx*w*Lw; 


V2=pi*Lc*(r4^2-r3 

-2); 

V3=pi*Lc*(r2^2-rl 

^2); 

V4=pi*ro^2*Lm; 



[m] Width of for and aft bulkhead. 

[m] Height of wall surrounding CCL module, 
[m] Height of inner/outer canister wall. 

[m] Length of missile. 

[m] Thickness of wall surrounding CCL module, 
[m] Near center of missile (1/16 in). 

[m] Radius of missile. 

[m] Inside radius of inner cylinder wall. 

[m] Outside radius of inner cylinder wall. 

[m] Inside radius of outer cylinder wall. 

[m] Outside radius of outer cylinder wall. 

[m] Distance between cylinder CL. 

[m/s^2] Gravitational constant. 

[W/m^2*KM] Stefan-Boltzmann constant. 
[m^2]Area of VLS wall perp to dir of heat flux. 

[m^2] Outer surface area of 1st cylinder. 

[m^2] Inner surface area of 1st cylinder. 

[m^2] Outer surface area of 2nd cylinder. 

[m^2] Inner surface area of 2nd cylinder. 

[m^2] Surface area of missile. 


% [K] INITIAL AMBIENT CONDITIONS TO START ITERATION. 
T2(l)=300; 
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T3(l)=300; 

T4(l)=300; 

T5(l)=300; 

T6(l)=300; 

T7(l)=300; 

T8(l)=300; 


% DEFINES TIME STEP AND NUMBER OF ITERATIONS 

S=60000; % Number of iterations. 

delt=0.36; % [s] Time step. 

% BEGINS LOOP TO SOLVE FOR ALL TEMPS IN RESISTANCE NETWORK. 

for p=l:S 
pl=p 


% [s] Time. 
t(p)=(delt*(p-l)); 

dT(p)=t(p)* 1.3333; 
dTA(p)=t(p)*0.0128; 
Tl(p)=300+dT(p); 
TlA(p)=1200+dTA(p); 

ifTl(p)<=1200 

TI(p)=TI(p); 

Tl(p)-TlA(p); 

end 

t(p+l)=(delt*(p)); 
dT(p+l)=t(p+l)* 1.3333; 
dTA(p+l)=t(p+l)*0.0128; 
Tl(p+l)=300+dT(p+l); 
TlA(p+l)=1200+dTA(p+l); 

ifTl(p+l)<=1200 

Tl(p+l)=Tl(p+l); 

else 

Tl(p+l)=TlA(p+l); 

end 
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%%% SOLVES SECTION 1 [T2] %%% 


% [K] Film temperature to evaluate properties of mild steel 
Tfl(p)=(Tl(p)+T2(p))/2; 

% [W/m*K] Thermal conductivity of bulkhead. 

kwl (p)=abs((l. 8869e-005 *Tfl (p)""2)-(6.9896e-002*Tfl (p))+83.43 8); 

% [J/kg*K] Specific Heat of bulkhead. 

cpl(p)=abs(((-2.2828e-006)*Tfl(p)^3)+((4.6812e-003)*Tfl(p)^2)-(2.3469*Tfl(p))+800.99); 

% [m^2/s] Thermal diffusivity of bulkhead. 

alfl(p)=abs(((4.0118e-014)*Tfl(p)^3)-((7.9337e-011)*Tfl(p)^2)+(3.3733e-008*Tfl(p))+1.0388 

e-005); 

% [kg/m'^3] Density of bulkhead, 
rho 1 (p)=kwl (p)/(cp 1 (p)*alfl (p)); 

% Solves for the conductive resistance in bulkhead. 

Rc 1 (p)=delx/(kw 1 (p) * A2); 

% [K] Film temperature to evaluate properties in 1st air space. 

Tf2(p)=(T2(p)+T3(p))/2; 

% [W/m*K] Thermal conductivity of air. 

ka2(p)=abs((4.9198e-011 *Tf2(p)^3)-(l .608e-007*Tf2(p)'^2)+(2.014e-004*Tf2(p))-0.020796); 

% [m^2/s] Kinematic viscosity of air. 

nu2(p)=abs(((7.786e-005*Tf2(p)^2)+(4.4053e-002*Tf2(p))-2.4972)*1.0e-006); 

% Prandtl #. 

Pr2(p)=abs((-2.5065e-011 *Tf2(p)'"3)+(8.4944e-008*Tf2(p)^2)-(l .0417e-004*Tf2(p))+0.74553); 

% [K^-1] Exp. Coefficient of air . 

B2(p)=l/Tf2(p); 

% Emissivity of mild steel. 

E2(p)=abs((-9.8416e-008*T2(p)''2)+(2.6912e-004*T2(p))+0.73699); 

% Emissivity of titanium. 

E3(p)=abs((-5.0412e-008*T3(p)^2)+(2.4949e-004*T3(p))+0.051504); 
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% Solves for Rayleigh and Grashof #'s WRT bulkhead inner surface. 
Gr2(p)=((g*B2(p)*(T2(p)-Tf2(p))*(Lw^3))/(nu2(p)^2)); 

Ra2(p)=Gr2(p)*Pr2(p); 

% Correlation to solve for the average Nusselt # over a 
% vertical flat plat. Recommended by Churchill and Chu, it is 
% good over entire range of Rayleigh #'s. 

Nu2(p)=(0.825+((0.387*(Ra2(p)^(l/6)))/((H-((0.492/Pr2(p))^(9/16)))^8/27))))^2; 

% Solves for convective resistance using the average convection 
% heat transfer coefficient over the bulkhead inner surface. 
h2(p)=Nu2(p)*ka2(p)/Lw; 

Rcola(p)=l/(h2(p)*A2); 

% Solves for Rayleigh and Grashof #'s WRT cylinder outer surface. 
Gr3(pHg*B2(p)*(T£2(p)-T3(p))*(Lc-3))/(nu2(p)^2); 

Ra3(p)=Gr3(p)*Pr2(p); 

% IF-THEN LOOP TO DETERMINE IF THE OUTER WALL OF THE CYLINDER 
% CAN BE APPROXIMATED AS A PLANE WALL. IF A PLANE WALL CHURCHILL 
% AND CHU'S CORRELATION FOR A VERTICAL SURFACE IS USED. IF 
% THERE ARE APPRECIABLE CURVATURE EFFECTS A CORRELATION DEVELOPED 
% BY LE FEVRE AND EDE BASED ON CYLINDER HEIGHT. 

ifGr3(p)=0 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)^2)/(5*(20+21*Pr2(p))))^0.25)+((4*Lc*(272+315*Pr2(p)))/(3 

5*2*r4*(64+63*Pr2(p)))); 

else 

if((2*r4)/Lc)>=35/(Gr3(p)^l/4) 

Nu3(p)=(0.825+((0.387*(Ra3(p)^(l/6)))/((l+((0.492/Pr2(p))^(9/16)))^(8/27))))^2; 

else 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)^2)/(5*(20+21*Pr2(p))))^.25)+((4*Lc*(272+315*Pr(p)))/(35 

*2*r4*(64+63*Pr2(p)))); 

end 

end 

% Solves for convective resistance using the average convection 
% heat transfer coefficient over outer surface of the cylinder. 
h3 (p)=(Nu3 (p)*ka2(p))/Lc; 

Rcolb(p)=l/(h3(p)*A3); 

% Adds convective resistances in series. 

Rco 1 (p)=Rco 1 a(p)+Rco 1 b(p); 
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% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, infinite plane to a row of pipes [C-6]. 
D=(2*r4)/s; 

F23=1 -((1 -D^2)^0.5)+(D*atan((( 1 -D^2)/D^2)^0.5)); 

% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

Rrl(pK(l-E2(p)V(E2(p)*A2))+(l/(A2*F23»+((l-E3(p)V(E3(p)*A3)); 

qol(p)=(sig*(T2(p)''4-T3(p)M))/Rrl(p); 


if qol(p)==0 

Rrad 1 (p)=Rr 1 (p)/(4 * sig * T3 (p)'^3); 
else 

Rradl(p)=l/((qol(p)/(T2(p)-T3(p)))+((4*sig*T3(pr3)/Rrl(^^^^ 

end 

% Solves for total resistance, [T3], and begins next iteration. 

Rt(p)=(Rco 1 (p)*Rrad 1 (p))/(Rrad 1 (p)+Rco 1 (p)); 

% Finite difference equation—explicit method. 

H1 (p)=delt/(rho 1 (p)*cp 1 (p)* V1); 

T2(p+l)=2*(((Hl(p)/Rcl(p))*(Tl(p)-T2(p)))+((Hl(p)/Rt(p))*(T3(p).T2(p)))+((T2(p)+Tl(p))/2) 

)-Tl(p+l); 

ifT2(p+l)<300 

T2(p+1)=300; 

else 

T2(p+l)=T2(p+l); 

end 


%%% SOLVES SECTION 2 [T3] %%% 

% [K] Film temperature to evaluate properties of titanium. 

Tf3(p)=(T3(p)+T4(p))/2; 

% [W/m*K] Thermal conductivity of canister wall. 

kc3(p)=abs((1.7011e-008*TD(p)^3)-(4.3316e-005*TD(p)^2)+(2.7279e-002*TD(p))+20.358); 
% [J/kg*K] Specific Heat of canister wall. 

cp3(p)=abs(((4.6480e-007)*TG(p)^3)-((1.4813e-003)*Tf3(p)^2)+(1.5387*Tf3(p))+139.93); 


73 





% [m^2/s] Thermal diffusivity of canister wall. 

alO(p)=abs(((-4.532e-015)*TD(p)'^3)+((l .7281e-01 l)*Tfi(p)^2)-(l .929e-008*Tf3(p))+l .3594e- 
005); 

% [kg/m^3] Density of canister wall. 
rho3(p)=kc3(p)/(alf3 (p)*cp3 (p)); 

% Solves for the conductive resistance through the cylinder wall. 
Rc2(p)=log(r4/r3)/(2*pi*kc3(p)*Lc); 

% Finite difference equation-explicit method. 

H2(p)=delt/(rho3(p)*cp3(p)*V2); 

T3(p+l)=2*(((H2(p)/Rt(p))*(T2(p)-T3(p)))+((H2(p)/Rc2(p))*(T4(p)-T3(p)))+((T3(p)+T4(p))/2) 

)-T4(p); 


%%% SOLVES SECTION 3 [T4] %%% 

% [K] Film temperature to evaluate air properties in 1st annulus. 

Tf4(p)=(T4(p)+T5(p))/2; 

% [W/m*K] Thermal conductivity of air. 

ka4(p)=abs((4.9198e-011 *Tf4(p)^3)-(l .608e-007*Tf4(p)^2)+(2.014e-004*Tf4(p))-0.020796); 

% [m^2/s] Kinematic viscosity of air. 

nu4(p)=abs(((7.786e-005*Tf4(p)^2)+(4.4053e-002*Tf4(p))-2.4972)*1.0e-006); 

% Prandtl #. 

Pr4(p)=abs((-2.5065e-011*Tf4(p)^3)+(8.4944e-008*Tf4(p)^2)-(1.0417e-004*Tf4(p))+0.74553); 

% [K^-1] Exp. Coefficient of air. 

B4(p)=l/Tf4(p); 

% Emissivity of titanium. 

E4(p)=abs((-5.0412e-008*T4(p)^2)+(2.4949e-004*T4(p))+0.051504); 

E5(p)=abs((-5.0412e-008*T5(p)^2)+(2.4949e-004*T5(p))+0.051504); 

% Solves for Rayleigh and Grashof #'s. 

Gr4(p)=((g*B4(p)*(T4(p)-T5(p))*(Lc^3))/(nu4(p)^2)); 

Ra4(p)=Gr4(p) *Pr4(p); 
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ifRa4(p)==0 

Ra4(p)=0.001; 

else 

Ra4(p)=Ra4(p); 

end 

% Nusselt correlation for curcular annuli heated and cooled on the 
% vertical curved surface. Developed by de Vahl Davis and Thomas. 
fpr4(p)=( 1 +((0.5/Pr4(p))^(9/l 6)))^(-16/9); 
Nu4(p)=0.364*((Ra4(p)*fpr4(p))^.25)*((r3/r2)^.5); 

% Solves for convective resistance using the average convection 
% heat transfer coefficient within the annuli. 
h4(p)=(Nu4(p)*ka4(p))/Lc; 

Rco2(p)=(A4+A5)/(h4(p)*A4*A5); 

% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, interior of outer right circular cylinder 
% of finite length to exterior of inner right circular coaxial 
% cylinder [C-87]. 

R=r3/r2; 

HH=Lc/r2; 

F45=(l/R)*(l-((HH'^2+R^2-l)/(4*HH))-l/pi*(acos((HH^2-R^2+l)/(HH^2+R^2-l))-(((sqrt(((H 

H^2+R^2+l)^2)-(2*(R^2))))/(2*HH))*acos((HH^2-R^2+l)/(R*(HH^2+R^2-l))))-(((HH^2-R^2 

+l)/(2*HH))*asin(l/R)))); 

% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

Rr2(p)=((l-E4(p))/(E4(p)*A4))+(l/(A4*F45))+((l-E5(p))/(E5(p)*A5)); 

qo2(p)=(sig*(T4(p)-4-T5(p)-4))/Rr2(p); 

if qo2(p)—0 

Rrad2(p)=Rr2(p)/(4*sig*T5(p)^3); 

else 

Rrad2(p)=l/((qo2(p)/(T4(p)-T5(p)))+((4*sig*T5(p)-3)/Rr2(p))); 

end 

% Solves for total resistance. 

Rt2(p)=((Rco2(p)*Rrad2(p))/(Rco2(p)+Rrad2(p))); 
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% Finite difference equation-explicit method. 

T4(p+l)=2*(((H2(p)/Rc2(p))*(T3(p).T4(p)))+((m(p)/Rt2(p))*(T5(p)-T4(p)))+((T3(p^^^ 

2)).T3(p+l); 


ifT4(p+l)<300 

T4(p+1)=300; 

0ls6 

T4(p+l)=T4(p+l); 

end 


%%% SOLVES SECTION 4 [T5] %%% 

% [K] Film temperature to evaluate properties of titanium. 

Tf5(p)=(T5(p)+T6(p))/2; 

% [W/m*K] Thermal conductivity of canister wall. 

lcc5(p)=abs((1.7011e-008*Tf5(p)^3)-(4.3316e-005*Tf5(p)^2)+(2.7279e-002*Tf5(p))+20.358); 

% [J/kg*K] Specific Heat of canister wall. 

cp5(p)=abs(((4.6480e-007)*Tf5(p)^3)-((1.4813e-003)*Tf5(p)^2)+(1.5387*Tf5(p))+139.93); 

% [m^2/s] Thermal diffiisivity of canister wall. 

alf5(p)=abs(((-4.532e-015)*Tf5(p)^3)+((1.728 le-01 l)*Tf5(p)^2)-(l .929e-008*Tf5(p))+l .3594e- 
005); 

% [kg/m'^3] Density of canister wall. 
rho5(p)=kc5(p)/(alf5(p)*cp5(p)); 

% Solves for the conductive resistance through the cylinder wall. 
Rc3(p)=log(r2/rl)/(2*pi*kc5(p)*Lc); 

% Finite difference equation—explicit method. 

H3(p)=delt/(rho5(p)*cp5(p)*V3); 

T5(p+l)=2*(((H3(p)/Rt2(p))*(T4(p)-T5(p)))+((H3(p)/Rc3(p))*(T6(p)-T5(p)))+((T5(p)+T6(p))/ 

2))-T6(p); 


%%% SOLVES SECTION 5 [T6] %%% 

% [K] Film temperature to evaluate air properties in 2nd annulus. 
TfB(p)=(T6(p)+T7(p))/2; 
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% [W/ni*K] Thermal conductivity. 

ka6(p)=abs((4.9198e-011*Tf6(p)^3)-(1.608e-007*Tf6(p)^2)+(2.014e-004*Tf6(p))-0.020796); 

% [m^2/s] Kinematic viscosity. 

nu6(p)=abs(((7.786e-005*Tf6(p)^2)+(4.4053e-002*Tf6(p))-2.4972)*1.0e-006); 

% Prandtl #. 

Pr6(p)=abs((-2.5065e-011*Tf6(p)^3)+(8.4944e-008*Tf6(p)'^2)-(1.0417e-004*Tf6(p))+0.74553); 

% [K^-1] Exp. Coefficient of air. 

B6(p)=l/Tf6(p); 

% Emissivity of titanium. 

E6(p)=abs((-5.0412e-008*T6(p)^2)+(2.4949e-004*T6(p))+0.051504); 

E7(p)=abs((-5.0412e-008*T7(p)^2)+(2.4949e-004*T7(p))+0.051504); 

% Solves for Rayleigh and Grashof #'s. 

Gr5(p)=(g*B6(p)*(T6(p)-T7(p))*(Lc^3))/(nu6(pr2); 

Ra5(p)=Gr5(p)*Pr6(p); 


ifRa5(p)==0 

Ra5(p)=0.001; 

else 

Ra5(p)=Ra5(p); 

end 

% Nusselt correlation for curcular annuli heated and cooled on the 
% vertical curved surface. Developed by de Vahl Davis and Thomas. 
fpr5(p)=( 1 +((0.5/Pr6(p))^(9/l 6)))^(-l 6/9); 
Nu5(p)=0.364*((Ra5(p)*fpr5(p))^.25)*((rl/ro)^0.5); 

% Solves for convective resistance using the average convection 
% heat transfer coefficient within the annuli. 
h5(p)=Nu5(p)*ka6(p)/Lc; 

Rco3(p)=(A6+A7)/(h5(p)*A6*A7); 

% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, interior of outer right circular cylinder 
% of finite length to exterior of inner right circular coaxial 
% cylinder [C-87]. 

Rl=rl/ro; 

HHl=Lc/ro; 
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F67=(l/Rl)*(l-((HH1^2+R1^2-l)/(4*HHl))-l/pi*(acos((HHl'^2-Rr2+l)/(HHr2+R1^2-l))-(( 

(sqrt(((HHK2+Rr2+l)^2)-(2*(R1^2))))/(2*HHl))*acos((HHr2-Rr2+l)/(Rl*(HHr2+Rr2 

-l))))-(((HH1^2-Rr2+l)/(2*HHl))*asin(l/Rl)))); 

% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

Rr3(p)=((l-E6(p))/(E6(p)*A6)>t(I/(A7*F67))+((l-E7(p)y(E7(p)*A7)); 

qo3(p)=(sig*Wpy'4-T7(p)''4))/Rr3{p); 

if qo3(p)==0 

Rrad3(p)=Rr3(p)/(4*sig*T7(p)'^3); 

else 

Rrad3(p)=l/((qo3(p)/(T6(p)-T7(p)))+((4*sig*T7(pr3)/Rr3(p))); 

end 

% Solves for total resistance. 

Rt3 (p)=(Rco3 (p)*Rrad3 (p))/(Rco3 (p)+Rrad3 (p)); 

% Finite difference equation—explicit method. 

T6(p+l)=2*(((H3(p)/Rc3(p))*(T5(p)-T6(p)))+((H3(p)/Rt3(p))*(T7(p)-T6(p)))+((T5(p)+T6(p))/ 

2))-T5(p+l); 

ifT6(p+l)<300 

T6(p+1)=300; 

else 

T6(p+l)=T6(p+l); 

end 

% [K] Film temperature to evaluate properties of aluminum. 

Tf7(p)=(T7(p)+T8(p))/2; 

% [W/m*K] Thermal conductivity of missile. 
kc7(p)=237; 

% [J/kg*K] Specific Heat of missile. 
cp7(p)=903; 

% [m^2/s] Thermal difiiisivity of missile. 
alf7(p)=97.1e-006; 

% [kg/m^3] Density of missile. 
rho7(p)=2702; 
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% Solves for the conductive resistance through the cylinder wall. 
Rc4(p)=log(ro/r)/(2*pi*kc7(p)*Lc); 

% Finite difference equation-explicit method. 

H4(p)=(delt/(rho7(p)*cp7(p)*V4)); 

T7(p+l)=2*(((H4(p)/Rt3(p))*(T6(p)-T7(p)))+((H4(p)/Rc4(p))*(T8(p)-T7(p)))+((T7(p)+T8(p))/ 

2))-T8(p); 

T8(p+l)=2*(((H4(p)/Rc4(p))*(T7(p)-T8(p)))+((T7(p)+T8(p))/2))-T7^ 

ifT8(p+l)<300 

T8(p+1)=300; 

else 

T8(p+l)=T8(p+l); 

end 

end 

t-t/60; 

plot(t,T8); grid; 
xlabel('Time (min)'); 
ylabel('Temperature (K)'); 
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APPENDIX C. CCL-FIRE PROGRAM CODE (2A-CORNER) 

This MATLAB code solves for the radiation resistance explicitly as described in 

the radiation section of Chapter II. This code is set up for scenario 2 and the comer CCL 

in which the adjoining bulkhead is assumed to be a reradiating surface. 

format bank 
clear 


% DIMENSIONS FOR CANISTERS & SURROUNDING VLS WALL. 


w=l; 

% 

wl=l; 

% 

Lw=l; 

% 

Lc=l; 

% 

Lm=l; 

% 

delx=0.0048; 

% 

r=0.00159; 

% 

ro=0.2586; 

% 

rl=0.2635; 

% 

r2=0.2667; 

% 

r3=0.3492; 

% 

r4=0.3571; 

% 

s=0.8382; 

% 

g=9.81; 

% 

sig=5.67e-008; 

% 

AR=wl*Lw; 

% 

A2=w*Lw; 

% 

A3=2*pi*r4*Lc; 

% 

A4=2*pi*r3*Lc; 

% 

A5=2*pi*r2*Lc; 

% 

A6=2*pi*rl*Lc; 

% 

A7=2*pi*ro*Lc; 

% 

V1 =delx* w*Lw; 


V2=pi *Lc*(r4^2-r3 

^2); 

V3=pi *Lc* (r2^2-r 1 

^2); 

V4=pi*ro^2*Lm; 



[m] Width of for and aft bulkhead. 

[m] Width of port and stbd bulkhead. 

[m] Height of wall surrounding CCL module. 

[m] Height of inner/outer canister wall. 

[m] Length of missile. 

[m] Thickness of wall surrounding CCL module, 
[m] Near center of missile (1/16 in). 

[m] Radius of missile. 

[m] Inside radius of inner cylinder wall. 

[m] Outside radius of inner cylinder wall. 

[m] Inside radius of outer cylinder wall. 

[m] Outside radius of outer cylinder wall. 

[m] Distance between cylinder CL. 

[m/s^2] Gravitational constant. 

[W/m^2*K^4] Stefan-Boltzman constant. 

[m^2] Area of reradiative bulkhead 

[m'^2] Area of VLS wall perp to dir of heat flux. 

[m^2] Outer surface area of 1st cylinder. 

[m^2] Inner surface area of 1st cylinder. 

[m^2] Outer surface area of 2nd cylinder. 

[m''2] Inner surface area of 2nd cylinder. 

[m^2] Surface area of missile. 


81 




% [K] INITIAL AMBIENT CONDITIONS TO START ITERATION. 


T2(I)=300; 

T3(l)=300; 

T4(l)=300; 

T5(l)=300; 

T6(1H00; 

T7(1^00; 

T8(l)=300; 

% DEFINES TIME STEP AND NUMBER OF ITERATIONS 

S=45000; % Number of iterations. 

delt=0.36; % [s] Time step. 


% BEGINS LOOP TO SOLVE FOR ALL TEMPS IN RESISTANCE NETWORK. 


for p=l:S 
pl=p 


% [s] Time. 
t(p)=(delt*(p-l)); 

% Temperature rise from ambient to 1325 K over a 5 minute period. 
% This temp is then maintained for the duration of the fire. 
dT(p)=t(p)*1.3333; 
dTA(p)=t(p)*0.0128; 

Tl(p)=300+dT(p); 

TlA(p)=1200+dTA(p); 

ifTl(p)<=1200 

Tl(p)=Tl(p); 

Tl(p)=TlA<p); 

end 

t(p+l)=(delt*(p)); 
dT(p+l)=t(p+l)* 1.3333; 
dTA(p+l)=t(p+l)*0.0128; 

Tl(p+l)=300+dT(p+l); 

TlA(p+l)=1200+dTA(p+l); 
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ifTl(p+l)<=1200 

Tl(p+l)=Tl(p+l); 

else 

Tl(p+l)=TlA(p+l); 

end 


%%% SOLVES SECTION 1 [T2] %%% 

% [K] Film temperature to evaluate properties of mild steel 
Tfl(pHTI(p)+T2(p))/2; 

% [W/m*K] Thermal conductivity of bulkhead. 
kwl(p)=abs((1.8869e-005*Tfl(p)^2)-(6.9896e-002*Tfl(p))+83.438); 

% [J/kg*K] Specific Heat of bulkhead. 

cpl(p)=abs(((-2.2828e-006)*Tfl(p)^3)+((4.6812e-003)*Tfl(p)^2)-(2.3469*Tfl(p))+800. 

99); 

% [m'^2/s] Thermal diffixsivity of bulkhead. 

alfl(p)=abs(((4.0118e-014)*Tfl(p)^3)-((7.9337e-011)*Tfl(p)^2)+(3.3733e-008*Tfl(p)) 

+1.0388e-005); 

% [kg/m^3] Density of bulkhead, 
rho 1 (p)=kwl (p)/(cp 1 (p)*alfl (p)); 

% Solves for the conductive resistance in bulkhead. 

Rc 1 (p)=delx/ (kw 1 (p) * A2); 

% [K] Film temperature to evaluate properties in 1st air space. 

Tf2(p)=(T2(p)+T3(p))/2; 

% [W/m*K] Thermal conductivity of air. 

ka2(p)=abs((4.9198e-011 *Tf2(p)^3)-(l .608e-007*Tf2(p)^2)+(2.014e-004*Tf2(p))-0.020 
796); 

% [m^2/s] Kinematic viscosity of air. 

nu2(p)=abs(((7.786e-005*Tf2(p)'^2)+(4.4053e-002*Tf2(p))-2.4972)*1.0e-006); 

% Prandtl #. 

Pr2(p)=abs((-2.5065e-011 *Tf2(p)^3)+(8.4944e-008*Tf2(p)^2)-(l .0417e-004*Tf2(p))+0. 
74553); 
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% [K^-1] Exp. Coefficient of air . 

B2(p)=l/Tf2(p); 

% Emissivity of mild steel. 

E2(p)=abs((-9.8416e-008*T2(p)^2)+(2.6912e-004*T2(p))+0.73699); 

% Emissivity of titanium. 

E3(p)=abs((-5.0412e-008*T3(p)^2)+(2.4949e-004*T3(p))+0.051504); 

% Solves for Rayleigh and Grashof #'s WRT bulkhead inner surface. 
Gr2(pH(g*B2(p)*(T2(p)-Tf2(p))*(Lw^3))/(nu2(pr^^^ 

Ra2(p)=Gr2(p)*Pr2(p); 

% Correlation to solve for the average Nusselt # over a 
% vertical flat plat. Recommended by Churchill and Chu, it is 
% good over entire range of Rayleigh #'s. 

Nu2(pH0.825+((0.387*(Ra2(p)^(l/6)))/((l+((0.492/Pr2(p))^(9/16)))^(8/27))))^2; 

% Solves for convective resistance using the average convection 
% heat transfer coefficient over the bulkhead inner surface. 
h2(p)=Nu2(p)*ka2(p)/Lw; 

Rco 1 a(p)= 1/(h2(p)* A2); 

% Solves for Rayleigh and Grashof #'s WRT cylinder outer surface. 
Gr3(p)=(g*B2(p)*(Tf2(p)-T3(p))*(Lc-3))/(nu2(p)-2); 

Ra3(p)=Gr3(p)*Pr2(p); 

% IF-THEN LOOP TO DETERMINE IF THE OUTER WALL OF THE CYLINDER 
% CAN BE APPROXIMATED AS A PLANE WALL. IF A PLANE WALL 
% CHURCHILL 

% AND CHU’S CORRELATION FOR A VERTICAL SURFACE IS USED. IF 
% THERE ARE APPRECIABLE CURVATURE EFFECTS A CORRELATION 
% DEVELOPED 

% BY LE FEVRE AND EDE BASED ON CYLINDER HEIGHT. 


ifGr3(p)=0 

Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)^2)/(5*(20+21*Pr2(p))))^0.25)+((4*Lc*(272+3I5*Pr2 

(p)))/(3 5 *2*r4*(64+63 *Pr2(p)))); 

else 

if ((2*r4)/Lc)>=3 5/(Gr3 (p^IM) 

Nu3(p)=(0.825+((0.387*(Ra3(p)^(l/6)))/((l+((0.492/Pr2(p))^(9/16)))^(8/27))))^2; 

else 
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Nu3(p)=((4/3)*((7*Gr3(p)*Pr2(p)^2)/(5*(20+21 *Pr2(p))))^0.25)+((4*Lc*(272+315*Pr2 

(p)))/(35*2*r4*(64+63*Pr2(p)))); 

end 

end 

% Solves for convective resistance using the average convection 
% heat transfer coefficient over outer surface of the cylinder. 
h3 (p)=(Nu3 (p)*ka2(p))/Lc; 

Rcolb(p)=l/(h3(p)*A3); 

% Adds convective resistances in series. 

Rco 1 (p)=Rco 1 a(p)+Rco 1 b(p); 

% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, infinite plane to a row of pipes [C-6]. 

D=(2*r4)/s; 

F23=l-((l-D^2)'^.5)+(D*atan(((l-D^2)/D^2)^.5)); 

FR3=F23; 

F2R=0.5 * (1 +(w 1 /w)-sqrt( 1 +((w 1 /w)^2))); 

% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

J1=1/(A2*F2R); 

J2=1/(AR*FR3); 

JT=(J1+J2)"^(-1); 

Rrl(p)=((l-E2(p))/(E2(p)*A2))+(l/((A2*F23)+JT))+((l-E3(p))/(E3(p)*A3)); 

Rradl(p)=l/((sigAlrl(p))*((T2(p)+T3(p))*(T2(p)-2+T3(pr^^^^ 

% Solves for total resistance, [T3], and begins next iteration. 

Rt(p)=(Rco 1 (p)*Rrad 1 (p))/(Rrad 1 (p)+Rco 1 (p)); 

% Finite difference equation—explicit method. 

H1 (p)=delt/(rho 1 (p) * cp 1 (p) * V1); 

T2(p+l)=2*(((Hl(p)/Rcl(p))*(Tl(p)-T2(p)))+((Hl(p)/Rt(p))*(T3(p)-T2(p)))+((T2(p)+T 

l(p))/2))-Tl(p+l); 

ifT2(p+l)<300 

T2(p+1)=300; 

else 

T2(p+l)=T2(p+l); 

end 


85 


%%% SOLVES SECTION 2 [T3] %%% 


% [K] Film temperature to evaluate properties of titanium. 

Tf3(p)=(T3(p)+T4(p))/2; 

% [W/m*K] Thermal conductivity of canister wall. 

kc3(p)=abs((1.7011e-008*TD(p)^3)-(4.3316e-005*Tf3(p)^2)+(2.7279e-002*Tfi(p))+20. 

358); 

% [J/kg*K] Specific Heat of canister wall. 

cp3(p)=abs(((4.6480e-007)*Tf3(p)^3)-((1.4813e-003)*Tf3(p)^2)+(1.5387*Tf3(p))+139.9 

3); 

% [m^2/s] Thermal diffusivity of canister wall. 

alf3(p)=abs(((-4.532e-015)*TD(p)^3)+((1.7281e-011)*Tf3(p)^2)-(1.929e-008*Tf3(p))+l 

.3594e-005); 

% [kg/m'^3] Density of canister wall. 
rho3(p)=kc3(p)/(alf3(p)*cp3(p)); 

% Solves for the conductive resistance through the cylinder wall. 
Rc2(p)=log(r4/r3)/(2*pi*kc3(p)*Lc); 

% Finite difference equation—explicit method. 

H2(p)=delt/(rho3(p)*cp3(p)*V2); 

T3(p+l)=2*(((H2(p)/Rt(p))*(T2(p)-T3(p)))+((H2(p)/Rc2(p))*(T4(p)-T3(p)))+((T3(p)+T 

4(p))/2))-T4(p); 


%%% SOLVES SECTION 3 [T4] %%% 

% [K] Film temperature to evaluate air properties in 1st annulus. 

Tf4(p)=(T4(p)+T5(p))/2; 

% [W/m*K] Thermal conductivity of air. 

ka4(p)=abs((4.9198e-011 *Tf4(p)^3)-(1.608e-007*Tf4(p)^2)+(2.014e-004*Tf4(p))-0 020 
796); 

% [m''2/s] Kinematic viscosity of air. 

nu4(p)=abs(((7.786e-005 *Tf4(p)^2)+(4.4053e-002*Tf4(p))-2.4972)* 1 .Oe-006); 
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% Prandtl #. 

Pr4(p)=abs((-2.5065e-011 *Tf4(p)^3)+(8.4944e-008*Tf4(p)^2)-(l .0417e-004*Tf4(p))+0, 
74553); 

% [K^-1] Exp. Coefficient of air. 

B4(p)=l/Tf4(p); 

% Emissivity of titanium. 

E4(p)=abs((-5.0412e-008*T4(p)^2)+(2.4949e-004*T4(p))+0.051504); 

E5(p)=abs((-5.0412e-008*T5(p)^2)+(2.4949e-004*T5(p))+0.051504); 

% Solves for Rayleigh and Grashof #'s. 
Gr4(pK(g*B4(p)*(T4(p).T5(p))‘(U-3)y(nu4(pr2)); 

Ra4(p)=Gr4(p)*Pr4(p); 

ifRa4(p)=0 

Ra4(p)=0.001; 

else 

Ra4(p)=Ra4(p); 

end 

% Nusselt correlation for circular annuli heated and cooled on the 
% vertical curved surface. Developed by de Vahl Davis and Thomas. 
fpr4(p)=(l+((0.5/Pr4(p))^(9/16)))^(-16/9); 
Nu4(p)=0.364*((Ra4(p)*fpr4(p))^0.25)*((r3/r2)^0.5); 

% Solves for convective resistance using the average convection 
% heat transfer coefficient within the annuli. 
h4(p)=(Nu4(p)*ka4(p))/Lc; 

Rco2(p)=(A4+A5)/(h4(p)*A4*A5); 

% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, interior of outer right circular cylinder 
% of finite length to exterior of inner right circular coaxial 
% cylinder [C-87]. 

R=r3/r2; 

HH=Lc/r2; 

F45=(l/R)*(l-((HH^2+R^2-l)/(4*HH))-l/pi*(acos((HH^2-R^2+l)/(HH^2+R^2-l))-(((sq 

rt(((HH'^2+R^2+l)^2)-(2*(R'^2))))/(2*HH))*acos((HH^2-R^2+l)/(R*(HH^2+R^2-l))))-( 

((HH^2-R^2+l)/(2*HH))*asin(l/R)))); 
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% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

Rr2(p)=((l-E4(p))/(E4(p)*A4))+(l/(A4*F45))+((l-E5(p))/(E5(p)*A5)); 

Rrad2(p)= l/((sig/Rr2(p))*((T4(p)+T5(p))*(T4(pr2+T5(pr2))); 

% Solves for total resistance. 

Rt2(p)=((Rco2(p)*Rrad2(p))/(Rco2(p)+Rrad2(p))); 

% Finite difference equation—explicit method. 

T4(p+l)=2*(((H2(p)/Rc2(p))*(T3(p)-T4(p)))+((H2(p)/Rt2(p))*(T5(p)-T4(p)))+((T3(p^^ 

T4(p))/2))-T3(p+l); 

ifT4(p+l)<300 

T4(p+1)=300; 

else 

T4(p+l)=T4(p+l); 

end 


%%% SOLVES SECTION 4 [T5] %%% 

% [K] Film temperature to evaluate properties of titanium. 

Tf5(pHT5(p)+T6(p)y2; 

% [W/m*K] Thermal conductivity of canister wall. 

kc5(p)=abs((l .701 le-008*Tf5(p)^3)-(4.3316e-005*Tf5(p)^2)+(2.7279e-002*Tf5(p))+20. 
358); 

% [J/kg*K] Specific Heat of canister wall. 

cp5(p)=abs(((4.6480e-007)*Tf5(p)^3)-((1.4813e-003)*Tf5(p)^2)+(1.5387*Tf5(p))+139.9 

3); 

% [m^2/s] Thermal difiusivity of canister wall. 

alf5(p)=abs(((-4.532e-015)*Tf5(p)^3)+((1.7281e-011)*Tf5(p)'^2)-(1.929e-008*Tf5(p))+l 

.3594e-005); 

% [kg/m^3] Density of canister wall. 
rho5(p)=kc5(p)/(alf5(p)*cp5(p)); 

% Solves for the conductive resistance through the cylinder wall. 
Rc3(p)=log(r2/rl)/(2*pi*kc5(p)*Lc); 
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% Finite difference equation—explicit method. 

H3(p)=delt/(rho5(p)*cp5(p)*V3); 

T5(p+l)=2*(((H3(p)/Rt2(p))*(T4(p)-T5(p)))+((H3(p)/Rc3(p))*(T6(p)-T5(p)))+((T5(p)+ 

T6(p))/2))-T6(p); 


%%% SOLVES SECTION 5 [T6] %%% 

% [K] Film temperature to evaluate air properties in 2nd annulus. 
Tf6(p)=(T6(p)+T7(p))/2; 

% [W/m*K] Thermal conductivity. 

ka6(p)=abs((4.9198e-011*Tf6(p)^3)-(1.608e-007*Tf6(p)^2)+(2.014e-004*Tf6(p))-0.020 

796); 

% [m'^2/s] Kinematic viscosity. 

nu6(p)=abs(((7.786e-005*Tf6(p)^2)+(4.4053e-002*Tf6(p))-2.4972)*1.0e-006); 

% Prandtl #. 

Pr6(p.)=abs((-2.5065e-011 *Tf6(p)^3)+(8.4944e-008*Tf6(p)^2)-(l .0417e-004*Tf6(p))+0. 
74553); 

% [K'^-1] Exp. Coefficient of air. 

B6(p)=l/Tf6(p); 

% Emissivity of titanium. 

E6(p)=abs((-5.0412e-008*T6(p)^2)+(2.4949e-004*T6(p))+0.051504); 

E7(p)==abs((-5.0412e-008*T7(p)^2)+(2.4949e-004*T7(p))+0.051504); 

% Solves for Rayleigh and Grashof #'s. 
Gr5(p)=(g*B6(p)*(T6(p)-T7(p))*(Lc-3))/(nu6(p)-2); 

Ra5(p)=Gr5(p)*Pr6(p); 

ifRa5(p)=0 

Ra5(p)=0.001; 

else 

Ra5(p)=Ra5(p); 

end 

% Nusselt correlation for circular annuli heated and cooled on the 
% vertical curved surface. Developed by de Vahl Davis and Thomas. 
fpr5(p)=(l+((0.5/Pr6(p))^(9/16)))^(-16/9); 
Nu5(p)=0.364*((Ra5(p)*fpr5(p))^0.25)*((rl/ro)^0.5); 
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% Solves for convective resistance using the average convection 
% heat transfer coefficient within the annuli. 
h5 (p)=Nu5 (p)*ka6(p)/Lc; 

Rco3(p)=(A6+A7)/(h5(p)*A6*A7); 

% Calculates view factor. Equation from the catalog of radiation 
% configuration factors, interior of outer right circular cylinder 
% of finite length to exterior of inner right circular coaxial 
% cylinder [C-87]. 

Rl=rl/ro; 

HHl=Lc/ro; 

F67=(l/Rl)*(l-((HH1^2+Rr2-l)/(4*HHl))-l/pi*(acos((HH1^2-Rr2+l)/(HHr2+Rl^ 

2-l)H((sqrt(((HH1^2+Rr2+l)^2)-(2*(RK2))))/(2*HHl))*acos((HHr2-R1^2+l)/(Rl* 

(HHr2+R1^2-l))))-(((HHr2-R1^2+l)/(2*HHl))*asin(l/Rl)))); 

% Solves the radiative resistance network for two surf enclosure, 

% and linearizes the resistance to be used ICW the convective 
% resistance. 

Rr3(p)=((l-E6(p))/(E6(p)*A6))+(l/(A7*F67))+((l-E7(p))/(E7(p)*A7)); 

Rrad3(p)= l/((sig/Rr3(p))*((T6(p)+T7(p))*(T6(pr2+T7(pr2))); 

% Solves for total resistance. 

Rt3 (p)=(Rco3 (p) *Rrad3 (p))/(Rco3 (p)+Rrad3 (p)); 

% Finite difference equation-explicit method. 

T6(p+l)=2*(((H3(p)/Rc3(p))*(T5(p)-T6(p)))+((H3(p)/Rt3(p))*(T7(p)-T6(p)))+((T5(p)+ 

T6(p))/2))-T5(p+l); 

ifT6(p+l)<300 

T6(p+1)=300; 

else 

T6(p+l)=T6(p+l); 

end 

% [K] Film temperature to evaluate properties of aluminum. 

Tf7(pHT7(p)+T8(p))/2; 

% [W/m*K] Thermal conductivity of missile. 
kc7(p)=237; 

% [J/kg*K] Specific Heat of missile. 
cp7(p)=903; 
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% [m^2/s] Thermal diffiisivity of missile. 
alf7(p)=97.1e-006; 


% [kg/m^3] Density of missile. 
rho7(p)=2702; 

% Solves for the conductive resistance through the cylinder wall. 

Rc4(p)=log(ro/r)/(2 * pi *kc7(p) *Lc); 

% Finite difference equation-explicit method. 

H4(p)=(delt/(rho7(p)*cp7(p)*V4)); 

T7(p+l)=2*(((H4(p)/Rt3(p))*(T6(p)-T7(p)))+((H4(p)/Rc4(p))*(T8(p)-T7(p)))+((T7(p)+ 

T8(p))/2))-T8(p); 

T8(p+l)=2*(((H4(p)/Rc4(p))*(T7(p)-T8(p)))+((T7(p)+T8(p))/2))-T7(^^ 

ifT8(p+l)<300 

T8(p+1)=300; 

else 

T8(p+l)=T8(p+l); 

end 

end 

t=t/60; 

plot(t,T8); grid; 
xlabel('Time (min)'); 
ylabel('Temperature (K)'); 
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APPENDIX D. TIME-TEMP PROFILE AT EACH NODE 

The following figures show the time-temperature profile for each node in the 
resistance network. These four figures reflect Scenario 1. Figures 24 and 25 reflect the 
CCL-FIRE program code which solves for surface radiative resistance explicitly. Figures 
26 and 27 reflect the CCL-FIRE program code which solves for the surface radiative 
resistance with a Taylor series expansion. 



Figure 24. Scenario 1A - Node Time-Temperature Profiles. 
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Figure 25. Scenario lA - Node Time-Temperature Profiles. 
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